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1.  Identification  and  Significance  of  the  Phase  I  project 

Modem  ship  designs  can  only  be  successful  if  the  ship  can  survive  in  a  hostile 
environment.  In  spite  of  susceptibility  measures  (measures  to  avoid  a  hit),  the 
importance  of  vulnerability  reduction  for  warships  cannot  be  emphasized  hard  enough.[l] 
Susceptibility  measures  can  never  eliminate  the  probability  of  a  hit  with  its  successive 
loss  of  lives,  casualties,  loss  of  mission  and  material.  External  threats  for  a  ship  become 
even  more  important  as  the  focal  point  for  naval  operations  has  shifted  towards  littoral 
areas,  where  ships  are  exposed  to  an  even  higher  risk. 

Engineers  are  required  to  make  assessments  and  create  new  ship  designs  that  will 
demonstrate  high  survivability.  Unfortunately,  current  analysis  methods  do  not  provide 
high  level  of  confidence,  and  instead  large  structural  design  safety  factors  are  employed. 
As  a  result,  ships  are  heavier  and  more  expensive  to  construct  and  maintain  than  may 
actually  be  required.  More  reliable  predictions  which  can  quantify  the  level  of 
confidence  and  can  incorporate  test  data  for  improving  the  confidence  in  the  predictions 
are  required.  Such  predictions  can  support  innovative  structural  configurations  in  new 
Navy  programs  such  as  the  DD(X),  and  provide  lower  life-cycle  costs  with  increased 
survivability. 

The  LS-DYNA  and  the  ABAQUS  commercial  codes  have  been  employed  for 
simulating  the  response  of  a  ship  structure  to  an  underwater  explosion.  The  LS- 
DYNA/USA  code  is  based  on  boundary  elements  and  simulates  the  underwater  shock 
fluid-structure  interaction  based  on  a  doubly  asymptotic  approximation.  It  models  the 
interaction  of  the  structure  with  the  surrounding  fluid  in  terms  of  wet  surface  variables 
and  it  eliminates  the  need  for  fluid  elements  around  the  outside  of  the  structure.  Coupled 
with  LS-DYNA  it  is  used  to  simulate  the  response  of  ships  and  submarines  to  shock  wave 
detonation  and  the  bubble  gas  effects.  DYSMAS  is  a  code  developed  by  NSWC/IHD 
specifically  for  modeling  accurately  the  damage  to  naval  structures  from  underwater 
explosions.  This  development  is  based  on  coupling  a  hydrocode  with  a  simulation  model 
of  the  target  structure  within  a  parallel  computing  environment. 

Up  to  date,  only  deterministic  simulations  have  been  performed  for  the  response  of 
ship  structures  subjected  to  an  explosive  detonation.  In  real  ship  structures  uncertainty 
exists  due  to  manufacturing  variability  and  imperfections  introduced  at  the  welds,  due  to 
variations  in  the  materials  (i.e.  properties  of  composites,  structural  damping  properties, 
stiffness  properties  due  to  residual  stresses  from  welding,  etc.).  Uncertainty  also 
originates  from  the  finite  element  and  boundary  element  modeling  practices  (i.e. 
modeling  the  connections  between  welded  plates,  including  all  the  scantlings  in  the 
model,  etc.)  and  from  capturing  accurately  the  load  from  the  shock  and  bubble  gas  effects 
exerted  on  the  ship  hull  (i.e.  accounting  properly  for  the  properties  of  the  mixed  phase 
medium  in  the  vicinity  of  the  ship  due  to  the  flow,  the  effects  from  the  characteristics  of 
the  sea-bed,  etc.).  The  results  from  a  deterministic  analysis  comprise  a  single  realization 
which  may  or  may  not  represent  accurately  the  actual  behavior  of  the  system  which  is 
been  modeled.  In  addition,  when  deterministic  models  are  employed  for  design 
refinement  and  optimization,  the  optimum  solution  typically  resides  at  the  boundaries  of 
the  feasible  region  with  respect  to  the  constraints.  Thus,  the  uncertainties  which  exist  in 
the  real  system  can  cause  the  optimized  solution  to  reside  in  the  infeasible  domain  and 
lead  into  a  degraded  and  undesirable  performance  for  the  optimum  design.  In  general, 
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performing  design  refinement  without  taking  into  account  uncertainty  can  be  misleading 
when  uncertainties  exist  in  the  real  physical  system  which  is  been  modeled. 

The  work  completed  during  the  Phase  I  effort  generated  a  new  capability  for  updating 
a  numerical  model  from  test  data  in  applications  of  ship  shock  analysis  due  to  underwater 
explosions,  and  for  providing  error  estimates  due  to  the  uncertainty  which  exists  in  the 
ship  structures,  the  modeling  techniques,  and  the  loads.  The  new  process  is  flexible 
enough  in  order  to  utilize  limited  measurements  which  may  be  available  for  model 
updating  and  it  considers  uncertainty  to  be  present  in  the  test  data.  The  end  result  is  an 
updated  model  with  a  higher  confidence  level  in  its  predictive  capabilities.  The 
confidence  level  associated  with  the  revised  simulation  model  is  also  part  of  the  final 
prediction.  The  new  development  works  with  non-linear  simulation  models  and  with 
any  variables  that  an  engineer  chooses  to  consider  as  random. 

The  Principal  Component  Analysis  (PC A)  [2-4],  the  Kriging  method  for  metamodel 
generation  [5-8],  concepts  of  validating  non-linear  finite  element  models  for  blast 
overpressure  inside  a  room  from  a  limited  set  of  test  data  [9],  and  principles  of  modeling 
uncertainty  and  predictive  accuracy  in  non-linear  finite  element  models  [10]  were 
employed  in  the  new  development. 

The  PCA  converts  the  time  histories  from  the  test  or  from  the  analysis  into  “modal” 
type  of  information.  The  term  “modal”  does  not  imply  that  there  are  any  linear 
limitations.  It  is  utilized  as  an  indication  that  the  time  domain  histories  are  decomposed 
into  more  fundamental  properties.  The  “modal”  information  is  used  for  constructing  fast 
running  models  which  are  employed  by  the  model  update  and  error  estimation 
procedures.  The  Kriging  method  is  used  to  generate  metamodels  for  the  principal 
components  of  the  numerical  models  and  for  constructing  the  fast  running  models.  The 
generation  of  the  fast  running  models  is  essential  in  performing  within  a  practical  amount 
of  time,  the  model  updating  under  uncertainty,  and  the  error  estimation,  using  large  finite 
element  models  from  commercial  codes.  The  Kriging  method  is  also  used  for  enriching 
the  test  data  if  only  a  limited  amount  of  data  is  available. 

The  following  list  of  tasks  was  completed  during  the  Phase  I  project.  The  related 
technical  information  is  presented  in  this  report: 

(i)  Metamodel  generation. 

(ii)  Principal  Component  Analysis  (PCA). 

(iii)  Enrichment  of  test  data. 

(iv)  Probabilistic  model  updating. 

(v)  Error  Estimation. 

(vi)  Case  studies: 

•  Section  of  generic  surface  ship 

•  Barge  for  which  test  data  and  a  numerical  model  were  provided  by 
Northrop  Grumman  Newport  News  (NN) 

During  the  Phase  I  effort  there  was  interaction  with  NN  and  with  NSWC/IHD.  The 
developments  of  the  Phase  I  project  established  a  solid  technical  foundation  for 
developments  and  validation  using  full  scale  ship  models.  Potential  topics  for  a  Phase  II 
effort  are  presented  at  the  end  of  this  report. 
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2.  Technical  Developments 
2.1  Metamodel  generation. 

A  metamodel  (or  surrogate  model)  is  comprised  by  equations  that  replace  the  actual 
simulation  model.  Metamodels  are  evaluated  from  actual  simulations  performed  at  a 
group  of  sample  points.  An  Optimal  Symmetric  Latin  Hypercube  (OSLH)  algorithm  that 
offers  a  combination  between  computational  effort  and  optimality  is  utilized  for 
identifying  the  sample  points.  Once  the  sample  points  and  the  values  for  the  performance 
variables  are  known,  the  Kriging  method  is  employed  for  generating  one  metamodel  for 
each  performance  function  of  interest.  The  Kriging  method  is  based  on  treating  Z  ( x ), 
the  error  between  the  actual  performance  variable  y  (5c )  and  an  approximate  value /? ,  as 
a  stochastic  process: 

y(x)=P  +  Z(x)  (1) 

where  x  is  the  d-dimensional  vector  of  the  random  variables  that  defines  the  point 
where  the  performance  variable  is  evaluated,  and  “d”  is  the  number  of  random  variables. 
Z(x  )  is  considered  as  a  normal  process  with  zero  mean  and  a  covariance  that  can  be 
expressed  as: 

cov  (Z  (x,),  Z  (*,.))  =  «t2R(x(,x;)  (2) 

where  a2  -is  the  process  variance  and  R(  xt ,  3cy )  is  the  spatial  correlation  function  (SCF). 
The  equation  used  for  the  SCF  is: 

R  (  x, ,  Xj )  =  n  exP(“  6  k  (x  a  -x,.*)2)  (3) 

*=1 

and  it  indicates  a  process  with  infinitely  differentiable  paths  in  the  mean  square  sense.  0  k 

is  the  correlation  parameter  that  corresponds  to  the  k"'  component  of  the  d-dimensional 
vector  of  the  random  variables  x ,  i.e.  k=l,2,..  .,d.  For  a  set  xs  comprised  by  “n”  number 

of  sampling  points  x  sj ,  i=l  ,2,. .  .,n 


xTs={x  sl,x  S2,...,X  sn}  (4) 

The  corresponding  performance  variable  ys  is  considered  known  and  its  values  are 
defined  as: 

y»  ={y  ( *  si )>  y  ( ■ *  s2 ),. •  •,  y( *  *)}  (5) 

The  vector  of  correlations  between  the  sample  points  xs  and  the  evaluation  point  x  can 
be  expressed  as: 

rr(r)={R(r,?J1),R(r,rj2),.,R(r,rffl)}  (6) 

The  correlation  matrix  [R]  is  also  defined  among  all  the  sample  points: 

[RMR  (*,,*,)]„  (7) 


The  spatial  correlation  function  in  Equations  (6)  and  (7)  has  been  defined  by  Equation 
(3).  In  the  Kriging  method  the  value  of  the  performance  function  evaluated  by  the 
metamodel  at  the  evaluation  point  x  is  treated  as  a  random  variable.  The  computation  of 
P  and  Z  ( x )  in  Equation  (1)  is  based  on  minimizing  the  mean  square  error  (MSE)  in  the 
response 
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MSE[j>(x)]=E[j>(*)-y(?)]2  (8) 

subjected  to  the  unbiasedness  constraint: 

E[jK*)]=E[y(x)]  (9) 

Thus,  a  set  of  parameters  0  k  (k=l  ,2,. .  .,d)  are  computed  as  the  solution  of  minimizing 

the  product  [(det(R))  •  <72  ], 

Wber e&2=  yn(ys-F  fi)T  R-'(ys-F  ft)  (10) 

F  is  a  vector  of  n-dimension  with  all  unit  entries  and  /?  defined  as: 

P  =  (FT  R-1  Fy1  FT  R-'ys  (11) 

Once  the  set  of  optimal  parameters  0  k  has  been  computed,  then  the  performance  function 
at  any  point  x  can  be  evaluated  as: 

y(x)  =  fi  +  rT(x)R~'(ys-F  $)  (12) 


where  rT  ( x )  is  defined  in  Equation  (6).  Equation  (12)  comprises  the  metamodel 
developed  by  the  Kriging  method. 

In  this  project  codes  based  on  the  OSLH  and  the  Kriging  methods  are  employed  for 
generating  all  of  the  metamodels.  The  computational  capability  of  generating 
metamodels  is  utilized  in  two  areas.  First,  metamodels  are  developed  for  enriching  the 
response  matrix  of  the  measurements  from  a  limited  amount  of  test  data.  Responses  can 
be  generated  for  points  on  the  ship  where  test  data  were  not  originally  collected. 
Metamodels  are  also  employed  during  the  probabilistic  model  updating  process  and  in  the 
assessment  of  the  error  estimation. 


2.2  Principal  Component  Analysis  (PCA) 

Principal  components  are  computed  for  the  response  matrix  from  the  test  and  from  the 
numerical  model.  Both  matrices  have  the  form: 


L  *«(*l)  -  Xm(0. 


(13) 


where  each  row  corresponds  to  a  different  location  either  on  the  model  or  on  the  physical 
system  where  responses  are  computed  or  measured,  respectively;  each  column 
corresponds  to  a  different  time  step.  The  PCA  allows  to  convert  the  rime  histories  into 
“modal”  type  of  information  through  a  singular  value  decomposition.  The  time  histories 
included  in  the  response  matrix  include  non-linear  effects,  such  as  yielding,  interaction 
between  the  load  and  the  flexible  structure,  etc.  The  response  matrix  [X]  is  written  as: 
m  =  [U]  [ W ]  [Vf  (14) 

where  [U]  is  a  matrix  which  contains  “modal”  type  of  information  about  the 
measurement  locations  in  each  column.  [W]  is  a  diagonal  matrix  which  contains 
information  about  the  energy  associated  with  each  mode  placed  in  each  column  of  [U]. 
[Vf  is  a  matrix  which  contains  time  domain  information  about  how  the  “modes”  are 
synthesized  in  order  to  produce  the  time  histories.  Figure  1  summarizes  the  information 
which  is  contained  in  the  principal  components. 
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m  Pfi 

'  '  Modal  participation 

Each  column  is  a  “mode"  Only  diagonal  terms  terms  at  each  time 

energy  in  each  “mode"  step 

Figure  1.  Information  contained  in  principal  components  of  a  response  matrix 

Principal  component  analysis  is  useful  for  non-linear  models  for  the  same  reason  that 
modal  properties  are  useful  in  linear  dynamic  systems  for  comparing  analysis  and  test 
data.  In  linear  structural  dynamics,  the  structural  modes  are  associated  with  resonant 
frequencies.  In  non-linear  systems  the  interpretation  of  modal  properties  in  general 
depends  on  the  selection  of  the  response  data  included  in  matrix  [X].  The  PC  A  provides 
a  compact  representation  of  the  response  for  a  non-linear  model.  Figure  2  presents  the 
first  four  “modes”  for  one  quarter  of  a  simply  supported  plate  with  an  impact  load  applied 
at  the  middle  of  the  plate.  It  can  be  observed  that  the  “modal”  information  contained  in 
each  one  of  the  columns  of  the  [U]  matrix  increases  in  complexity  for  higher  modes  and 
that  each  “mode”  satisfies  the  physical  boundary  conditions. 


Figure  2.  Example  of  PCA  “modes”  for  the  quarter  of  a  simply  supported  plate  due  to  an 
impact  load  applied  at  the  middle  of  the  plate 
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In  order  to  perform  the  principal  component  decomposition,  equation  (14)  is  re-written 
as: 


m=[®  <t>  j 

where  [  W ]  = 


~D 

o' 

V 

0 

0 

fir  _ 

)  O' 

I  0 

,[U]~- 

*,] ,  w?  = 


V 

Vr 


(15) 


The  principal  components  decomposition  is  based  on  equation  (15)  and  results  in  the 
following  expression  for  [X]: 

[X]  =  m[D][?]]  (16) 

Equation  (16)  constitutes  the  principal  components  decomposition  of  the  response  matrix 
either  from  the  test  or  from  the  analysis.  The  matrices  [<E>J  and  [?]T]  do  not  participate 

in  the  definition  of  matrix  [X]  since  they  are  both  multiplied  with  the  zero  entries  of 
matrix  [W],  and  therefore  they  are  not  considered  in  the  definition  of  the  principal 
components. 

In  this  project  a  code  has  been  developed  for  performing  PCA  analysis  for  a  response 
matrix.  The  PCA  code  is  employed  for  probabilistic  model  update  and  for  error 
estimation.  The  corresponding  results  demonstrate  the  utilization  of  the  PCA  code. 


2.3  Enrichment  of  test  data 

Since  test  data  only  at  a  limited  number  of  locations  may  be  available,  a  capability  to 
essentially  increase  the  number  of  rows  in  the  response  [X]  matrix  was  developed.  For 
each  time  step  the  physical  information  in  a  column  of  the  test  matrix  is  used  along  with 
the  corresponding  coordinates  of  the  locations  where  the  test  data  were  collected,  in  order 
to  generate  metamodels  for  the  measured  data.  The  Kriging  code  is  employed  here  for 
enriching  the  test  data.  The  process  is  summarized  in  Figure  3 : 

Enriched  test  matrix 


Original  test  matrix 

*i(0 

•  *,(0 

XO  -  *i(0" 

XM  • 

•  XO 

3.(0  -  XO 

Enrichment  of  data 
through  metamodels 

j^wi+m'CO 

Figure  3.  Process  of  enriching  the  test  data 


Results  form  two  different  applications  are  presented  in  order  to  validate  the  capability  of 
enriching  the  test  data.  The  first  application  is  associated  with  a  rectangular  plate 
subjected  at  its  center  to  an  impact  load  of  short  duration.  The  finite  element  model  of 
this  plate  has  -2,700  degrees  of  freedom  (dof)  and  due  to  its  relative  small  size,  this 
model  was  utilized  for  the  initial  validation  of  all  new  developments.  The  duration  of  the 
impact  load  is  selected  in  such  a  manner  to  represent  typical  durations  from  shock  loads 
in  naval  applications.  The  finite  element  model  for  the  plate  is  depicted  in  Figure  4  and 
the  load  time  history  is  depicted  in  Figure  5.  LS-DYNA  is  utilized  for  the  analysis. 
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Figure  4.  Finite  element  model  for  the  plate  utilized  in  the  validation/demonstration. 


s 


Figure  5.  Time  history  of  the  load  applied  on  the  center  of  the  plate. 

The  second  application  is  associated  with  a  section  of  a  typical  surface  ship  subjected  to  a 
load  with  time  history  similar  to  typical  loads  in  underwater  explosion  applications.  The 
finite  element  model  for  a  typical  ship  is  presented  in  Figure  6.  LS-DYNA  is  employed 
again  for  performing  the  non-liner  FEA  analyses.  Each  finite  element  in  the  ship  model 
has  approximate  main  dimension  of  1ft.  The  ship  section  is  analyzed  as  part  of  the  case 
study  of  this  project.  All  the  steps  which  were  originally  validated  and  demonstrated 
through  the  simple  plate  model  are  validated  and  demonstrated  as  well  by  the  ship  section 
model  in  a  parallel  effort.  The  ship  section  model  and  the  corresponding  loads  were 
selected  in  collaboration  with  NN  in  order  to  make  certain  that  the  case  study  has  Naval 
relevance.  The  ship  section  finite  element  model  has  ~5 5,000  dof.  The  model  contains 
bottom  frame  members,  stiffeners,  bulkheads,  and  a  double  bottom  structure. 
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Figure  6.  Typical  ship  model  (left),  and  the  section  of  the  ship  model  (right)  which  is 

utilized  in  the  case  study 

In  order  to  validate  and  demonstrate  how  the  enrichment  of  the  test  data  process  operates 
15  points  marked  in  Figure  4  with  dark  circles  are  used  as  sample  points  for  developing 
metamodels  for  the  response  at  the  6  points  which  are  marked  with  plain  circles  in  Figure 
4.  Analysis  is  performed  by  LS-DYNA.  The  time  histories  at  the  15  sample  points 
comprise  the  original  response  matrix.  The  Kriging  code  is  used  to  generate  the  time 
histories  for  the  response  at  the  other  6  points  and  augment  the  original  response  matrix, 
thus  generating  the  enriched  response  matrix.  The  time  histories  computed  by  the 
Kriging  code  are  compared  to  the  actual  LS-DYNA  results  (Figure  7).  Very  good 
correlation  is  observed  between  the  actual  LS-DYNA  results  and  the  acceleration  time 
histories  created  by  the  Kriging  code  for  all  6  evaluation  points.  The  correlation  for  two 
representative  nodes  is  presented  in  Figure  7. 


Figure  7.  Time  histories  for  the  enriched  response  matrix  computed  by  the  Kriging  code 

and  compared  to  actual  LS-DYNA  results 

A  similar  analysis  is  performed  using  the  ship  section  (Figure  6).  Since  the  induced 
velocity  response  at  bulkheads  is  a  primary  parameter  measured  in  a  typical  underwater 
explosion  test,  the  normal  velocity  at  the  12  nodes  marked  with  dark  circles  in  Figure  6 
comprise  the  information  at  the  sample  points.  Time  histories  for  the  induced  velocity  at 
the  six  points  marked  with  the  plain  circles  in  Figure  6  are  generated  by  the  Kriging  code 
from  the  information  at  the  sample  points.  The  time  histories  generated  by  the  Kriging 
code  and  those  computed  by  the  actual  LS-DYNA  analysis  are  compared  for  all  6 
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evaluation  points.  Results  for  3  evaluation  points  (one  from  each  bulkhead)  are  presented 
in  Figure  8. 


Node-172  Node-806 


Time  (s) 

Figure  8.  Time  histories  for  the  enriched  response  matrix  computed  by  the  Kriging  code 
and  compared  to  actual  LS-DYNA  results  for  the  case  study  model 

Very  good  correlation  is  observed  for  the  velocities  at  all  the  points  between  the  Kriging 
code  and  the  actual  LS-DYNA  analysis.  As  a  clarification,  the  LS-DYNA  results  at  the 
evaluation  points  were  not  utilized  when  computing  the  response  at  the  evaluation  points 
through  the  Kriging  code.  Only  the  LS-DYNA  results  at  the  sample  points  were 
employed  for  this  purpose  in  both  applications. 

In  order  to  further  demonstrate  and  validate  the  enrichment  process,  PCA  is  performed 
for  a  response  matrix  comprised  exclusively  from  LS-DYNA  results  and  a  response 
matrix  comprised  by  LS_DYNA  results  for  some  points  and  by  enriched  data  for  other 
points.  The  plate  presented  in  Figure  9  is  utilized  again 
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Figure  9.  Sample  points  and  evaluation  points  for  used  in  PCA  and  Kriging  code 

validation 

In  this  case  the  LS-DYNA  results  at  the  9  points  marked  with  the  dark  circles  along  with 
the  time  histories  computed  by  the  Kriging  method  at  the  6  points  marked  with  the  plain 
circles  are  used  for  defining  the  response  matrix.  The  PCA  is  performed  for  the  response 
matrix.  The  PCA  results  are  compared  with  the  PCA  results  computed  for  a  response 
matrix  which  is  comprised  by  all  1 5  (9+6)  time  histories  computed  by  LS-DYNA.  Since 
the  same  PCA  results  are  generated  for  both  response  matrices,  this  indicates  that  the 
enrichment  of  the  data  is  performed  correctly.  Using  the  PCA  results,  the  vibration  at  the 
6  points  marked  with  the  plain  circles  is  also  reconstructed  and  compared  with  the  actual 
LS-DYNA  vibration  responses  at  these  points  (Figure  10). 


Figure  10.  Reconstructed  vibration  from  PCA  when  using  the  time  histories  from  both 
LS-DYNA  and  from  the  Kriging  method  vs.  the  actual  LS-DYNA  vibration  results 

The  good  correlation  between  the  reconstructed  responses  in  Figure  10,  demonstrates  that 
the  Kriging  code  provides  good  quality  time  histories  and  that  the  concept  of  enriching 
the  test  data  response  matrix  works  properly. 

Similar  PCA  and  vibration  reconstruction  from  PCA  results  is  performed  for  the  ship 
section  of  the  case  study  (Figure  6).  The  LS-DYNA  response  at  the  12  points  marked 
with  dark  circles  and  the  Kriging  results  at  the  6  points  marked  with  plain  circles  in 
Figure  6  constitute  one  response  matrix,  and  the  corresponding  results  are  labeled  as 
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“Kriging”  in  the  Figures.  The  LS-DYNA  results  at  all  18  (12+6)  points  constitute  a 
second  response  matrix.  PCA  are  performed  for  both.  The  diagonal  terms  of  matrix  [D] 
(Equation  16),  representing  the  energy  in  each  PCA  “mode”,  are  presented  in  Figure  1 1 
for  both  analyses. 


Figure  11.  Diagonal  terms  of  [D]  matrix  from  PCA  analyses 


The  reconstructed  response  from  the  principal  components  created  by  both  PCA  analyses 
are  presented  in  Figure  12  for  two  of  the  nodes. 


8PCs-Node:172  8PCs-Node:806 


Time  (s)  Time  (s) 


Figure  12.  Reconstructed  response  from  the  PCA  analyses 


Very  good  correlation  is  observed  between  all  the  results  produced  from  the  PCA 
analysis  based  on  all  LS-DYNA  results,  and  the  PCA  analysis  which  is  based  on  some 
LS-DYNA  results  and  some  Kriging  results.  Thus,  once  again  the  test  data  enrichment 
process  and  the  PCA  development  are  validated  through  the  case  study  results. 

2.4  Probabilistic  model  updating 

Model-to-test  correlation  metrics  similar  to  the  ones  used  in  linear  structural  dynamics 
are  established  based  on  the  principal  components  of  the  test  and  the  response  matrices. 
The  following  set  of  model-to-test  correlation  metrics  are  defined: 

A'P  =  0-7  (17  a) 
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the  left  superscript  zero  indicates  variables  from  the  analysis  and  variables  without 
superscript  originate  from  the  principal  components  of  the  test  matrix.  The  measures 
established  by  equations  (17a  -  17c)  are  generic.  The  cross  orthogonality  matrices  on  the 
right  hand  side  of  equations  (17a)  and  (17b)  are  the  directional  cosine  matrices  between 
the  respective  spaces  spanned  by  the  principal  vectors.  The  scaled  differences  of  singular 
values  in  equation  (17c)  represent  the  differences  in  energy  content  for  the  corresponding 
principal  components.  The  uncertainty  of  the  system  response  is  determined  by 
considering  the  variation  of  the  response  as  a  function  of  the  normalized  model-to-test 
correlation  matrices. 

The  variation  in  the  correlation  metrics  between  the  test  and  analytical  results  is 
defined  as: 


A r 


'  vec(A'F) ' 

•  vec(Av)  > 
diag(AD) 


(18) 


where  the  symbols  vec( )  and  diag( )  in  equation  (18)  represent  respectively,  a  column¬ 
wise  vectorization  of  a  matrix  and  a  vector  formed  by  the  diagonal  terms  of  a  matrix. 

By  considering  uncertainty  in  the  parameters  y  of  the  numerical  model,  equation  (16) 
becomes: 

(i») 

where  the  columns  of  [4>(y)]  and  the  diagonal  entries  of  [£>(y)]  are  evaluated  from 
metamodels  developed  by  the  Kriging  method.  Equation  (19)  represents  a  fast-running 
model  for  the  response  predicted  by  the  analytical  model.  Utilizing  fast-running  models 
is  critical  in  performing  the  model  updating  and  the  error  estimation  for  any  practical 
finite  element  model. 

A  common  approach  for  minimizing  the  error  between  model  predictions  and 
experimental  data,  while  introducing  statistical  information  into  the  process  is  through 
minimization  of  a  Bayesian-type  of  objective  function: 

j  =  (x-  [£„]-'  (l  -  °x(r))T  +  (r„  -  r)  [s„]  (r,  -  rf  m 


where  [S1^]  is  the  covariance  matrix  of  measurement  uncertainty,  [S^,]  is  the 
covariance  matrix  of  parametric  uncertainty  in  the  model,  yo  is  the  vector  of  the  original 


parameters,  y  is  the  updated  vector  of  parameters,  X  is  the  vectorized  measured  time 
history,  and  °X{y)  is  the  vectorized  time  history  from  the  simulation  as  it  is  evaluated 
from  equation  (19).  The  covariance  matrices  [S^]  and  [5^]  embody  the  a  priori 

knowledge  and  beliefs  regarding  the  errors  in  the  measurements  and  the  uncertainty  in 
the  parameters  of  the  numerical  model.  The  corresponding  Fisher  information  matrices, 
[Fja]  and  [F^,]  are  the  pseudo-inverse  of  [5^ ]  and  [S^],  respectively.  The 

covariance  matrices  act  as  weighting  factors  placing  more  importance  on  parameters 
associated  with  higher  levels  of  certainty.  The  first  term  in  the  functional  defined  by 
Equation  (20)  minimizes  the  difference  between  the  test  data  and  the  numerical  results. 
The  second  term  in  Equation  (20)  imposes  a  requirement  for  a  minimal  possible  change 
from  the  original  numerical  model  during  the  model  update  process.  The  latter 
requirement  contributes  to  the  uniqueness  of  the  solution  from  the  optimization  process. 
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Once  the  new  parameter  estimates  are  derived,  the  parameter  covariance  matrix  [5^] 
is  updated  through  the  first  order  Bayesian  method: 


fc]= 

where 


KY + [tJ  \p„y  [rXrf 


(21) 


-S^]  is  the  revised  covariance  matrix  of  the  parameter  error,  and  \rxy\  is  the 
sensitivity  matrix.  The  updated  pseudo-inverse  matrix  [f^]  is  evaluated  as: 

KbkJ+KT  farlfcJ  (22) 


An  optimization  code  has  been  developed  based  on  Equations  (20),  (21),  and  (22)  and  is 
used  for  updating  the  simulation  model.  This  updating  process  takes  into  account  the 
confidence  in  the  test  data  and  the  uncertainties  in  the  numerical  model.  Thus,  the 
updated  numerical  model  correlates  the  best  to  the  test  data  in  a  probabilistic  sense.  It  is 
essential  to  update  the  model  first  in  a  probabilistic  sense  before  employing  it  for 
analyses  which  include  error  estimation  in  the  prediction  in  order  to  improve  the 
confidence  in  the  numerical  results. 

In  order  to  validate  and  demonstrate  the  model  update  process  the  plate  model 
presented  in  Figure  9  is  utilized.  Concentrated  masses  equaling  the  35%  of  the  total  mass 
of  the  plate  are  distributed  over  the  plate.  The  magnitude  of  the  load,  the  yielding  stress, 
and  the  values  of  the  distributed  masses  were  perturbed  from  the  configuration  that 
represents  the  “test”  data.  Table  1  summarizes  the  values  of  the  variables  assigned  to  the 
configuration  representing  the  “test”  data,  the  original  numerical  model,  mid  the  standard 
deviation  assigned  to  each  one  of  the  variables  in  order  to  incorporate  the  uncertainty  of 
the  numerical  model  in  the  computations. 


“Test” 

Original  Model 

Sigma 

Yielding  strength 

hesbsbi 

■EBB 

951 

Magnitude  of  load 

5.0E4  N 

4.0E4N 

0.5E4N 

Table  1.  Summary  of  values  for  variables  with  uncertainty  in  pi 

ate  impact  model 

For  the  three  points  identified  in  Figure  13,  the  “test”  data,  the  results  corresponding  to 
starting  point  for  the  model  update  process,  and  the  results  corresponding  to  the  update 
model  are  presented  in  Figure  14.  As  it  can  be  observed,  the  perturbation  which  was 
introduced  in  the  model  resulted  in  a  large  initial  difference  between  the  “test”  data  and 
the  starting  point  of  the  model  update  process.  However,  after  the  model  has  been 
updated  the  “test”  data  and  the  results  for  the  updated  model  correlate  well. 
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Figure  13.  Three  points  for  which  results  are  presented  before  and  after  the  model  update 
process 


Figure  14.  Response  time  histories  for  the  “test”  data,  the  model  used  as  starting  point  in 
the  model  updating  process  and  the  finally  updated  model 

Initially  only  the  fast  running  model  (Equation  19)  was  utilized  during  the  model  update 
process.  Although  the  starting  point  for  the  model  update  process  differed  significantly 
from  the  “test”  data,  the  updated  model  correlated  well  with  the  “test”  data.  This 
demonstrates  that  the  probabilistic  model  updating  process  can  handle  large  initial 
differences  between  the  numerical  model  and  the  test  data,  while  accounting  for  the 
uncertainty  in  the  test  and  in  the  numerical  model.  A  different  scale  is  used  for  each  one 
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of  the  three  plots  in  Figure  14  since  points  closest  to  the  impact  location  exhibit  higher 
response. 

The  model  update  procedure  was  performed  for  a  second  time.  The  actual  LS-DYNA 
code  was  executed  every  time  that  the  model  updating  program  required  a  new  evaluation 
of  the  structural  response.  The  final  results  for  the  updated  model  also  matched  the  “test” 
data  during  the  second  model  update  analysis.  However,  a  significant  time  reduction  was 
achieved  when  using  the  fast  running  model  as  it  can  be  observed  from  the  cpu  time 
comparison  presented  in  Table  2. 


Run  Time 

LS-DYNA 

~1 1  hours 

Fast  running  model 

~12  minutes 

Table  2.  cpu  time  comparison  between  model  update  based  on  LS-DYNA  and  model 

update  based  on  a  fast  running  model 

The  results  demonstrate  that  the  model  update  code  operates  correctly,  that  the  Kriging 
code  works  properly  when  generating  the  fast  running  model,  and  that  significant  time 
savings  can  be  realized  when  using  the  fast  running  model.  Utilizing  a  fast  running 
model  in  the  model  update  process  is  crucial  since  the  large  LS-DYNA  models  used  in 
ship  applications  would  otherwise  make  the  model  update  process  infeasible. 

2.5  Error  Estimation 

Once  the  numerical  model  has  been  updated,  the  fast  running  model  based  on 
Equation  19  and  the  principal  component  matrices  [D(y)]  and  [r|(Y)]  computed  at  the  state 
of  the  updated  model  are  employed  for  assessing  the  standard  deviation  in  the  response. 
Due  to  the  small  cpu  time  requirements  of  the  fast  running  models  a  Monte  Carlo 
simulation  is  performed  in  order  to  generate  information  about  the  standard  deviation  of 
the  response  from  the  uncertainty  in  the  parameters  of  the  model.  The  expected  standard 
deviation  in  the  response  is  computed  for  the  updated  model  of  the  plate  presented  in  the 
previous  Section  (Figure  14).  Results  for  the  velocity  time  histories  from  the  “test”  data, 
the  original  numerical  model,  the  updated  numerical  model,  and  +  or  -  two  standard 
deviations  from  the  updated  model  are  presented  in  Figure  15.  In  this  error  estimation 
computations  the  standard  deviation  in  the  parameters  with  uncertainty  was  obtained 
from  Table  1.  The  standard  deviation  of  the  response  depends  on  the  uncertainty 
considered  in  the  parameters  of  the  numerical  model. 
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Figure  15.  Response  time  histories  for  the  “test’  data,  the  starting  model  in  the  model 
update  process,  the  updated  model,  and  the  +  or  -  two  standard  deviations 

Due  to  the  intentionally  large  difference  between  the  starting  point  for  the  model  update 
process  and  the  “test’  data,  the  line  which  corresponds  to  the  original  model  is  often 
beyond  the  +  or  -  two  standard  deviations.  However,  the  difference  between  the  updated 
model  and  the  “test”  data  is  almost  exclusively  within  the  +  or  -  two  standard  deviations. 
The  total  run  time  for  the  error  estimation  was  of  the  order  of  1  cpu  min  on  a  Pentium  4 
PC. 

2.6  Case  study 

It  is  important  to  demonstrate  that  the  new  developments  can  handle  the  probabilistic 
model  update  and  the  error  estimation  with  finite  element  models  that  contain  a  large 
number  of  dof,  and  that  they  can  interact  with  commercial  finite  element  codes.  These 
capabilities  are  demonstrated  through  the  case  studies.  The  probabilistic  model  update 
and  the  error  estimation  capabilities  presented  in  Sections  2.4  and  2.5  were  applied  on 
two  different  models.  First,  The  model  of  the  ship  section  presented  in  Figure  6  is  used 
for  validating  and  showcasing  the  developed  capabilities  on  a  model  significantly  larger 
and  more  demanding  than  the  simple  plate  model  (Figure  13)  which  was  utilized  during 
the  development  stages  of  the  model  update  and  error  estimation  process.  A  real 
application  from  an  UNDEX  test  of  a  barge  and  the  corresponding  ABAQUS  model  were 
provided  by  NN  and  utilized  to  validate  the  new  developments  of  this  SBIR  project.  This 
barge  application  was  not  part  of  the  original  scope  of  the  project.  However,  an  extra 
effort  was  made  by  MES  in  order  to  demonstrate  that  the  new  developments  work  for  real 
applications  of  Naval  interests. 
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2.6.1.  Probabilistic  model  update  and  error  estimation  for  a  hull  section  of  a  generic 
ship 

The  section  of  the  ship  hull  used  in  this  analysis  is  depicted  in  Figure  6.  The 
velocities  at  the  points  on  the  longitudinal  bulkheads  which  are  identified  in  Figure  6  are 
considered  in  the  model  update  process.  The  distribution  of  concentrated  masses 
supported  by  the  deck  structure,  and  the  magnitude  of  the  impact  load  applied  at  the 
bottom  structure  are  considered  as  the  parameters  with  uncertainty.  Table  3  summarizes 
the  values  that  these  parameters  acquire  for  the  “test”  configuration,  for  the  starting  point 
of  the  model  update  process,  and  for  the  standard  deviation  considered  in  the  variables. 


Loading  (N) 

Mass  (Kg) 

7.5E6 

3750 

1.5E7 

4125 

Standard 

deviation 

1.125E6 

960.5 

Table  3.  Summary  of  values  for  variables  with  uncertainty  in  the  ship  hull  analysis 

The  value  of  the  load  was  selected  high  enough  in  order  to  introduce  yielding  in  the 
longitudinal  bulkheads.  Figure  16  presents  results  for  three  points,  one  on  each 
longitudinal  bulkhead  (Figure  6).  For  each  point  the  results  from  the  “test”  data,  from  the 
configuration  before  the  model  update,  from  the  updated  model,  and  from  +  or  -  two 
standard  deviations  away  from  the  updated  model  are  presented. 


node-169  node-812 


Time  (s) 

Figure  16.  Initial  difference  between  “test”  (red  line)  and  numerical  model  (blue  line)  for 

configuration  used  in  the  case  study 
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As  it  can  be  observed  the  results  from  the  original  model  (which  comprises  the  starting 
point  for  the  model  updating  process)  differ  significantly  from  the  “test”  results.  This 
was  done  intentionally  in  order  to  demonstrate  that  the  model  update  process  works 
properly  in  adjusting  the  model  in  order  to  match  the  test  data.  The  updated  model 
matches  the  “test”  data  well,  and  the  difference  between  the  updated  model  and  the  test 
data  is  smaller  than  the  two  standard  deviations  margin  which  is  also  plotted  in  the 
Figures.  In  order  to  place  the  cpu  requirements  of  the  model  update  and  of  the  error 
estimation  process  into  prospective,  a  single  complete  LS-DYNA  simulation  for  the 
section  of  the  ship  hull  requires  ~3  cpu  hours.  Once  the  fast  running  model  has  been 
generated,  the  probabilistic  model  update  process  requires  ~25  cpu  min,  and  the  error 
estimation  process  requires  ~6  cpu  min  on  a  Pentium  III  PC.  Thus,  the  fast  running 
models  which  are  based  on  PCA  make  these  analyses  feasible  for  large  finite  element 
models. 

2.6.2.  Probabilistic  model  update  and  error  estimation  for  a  barge 

The  model  of  a  barge  in  ABAQUS  was  provided  by  Jay  Warren  of  NN.  A  picture  of 
the  barge  is  presented  in  Figure  17  [1 1].  The  ABAQUS  model  contained  both  structural 
elements  and  fluid  elements.  The  finite  element  model  contains  4,934  structural  elements 
and  74,560  fluid  elements.  The  complete  finite  element  model  is  depicted  in  Figure  18, 
and  the  structural  portion  of  the  finite  element  model  in  Figure  19.  The  charge  was 
prescribed  as  a  high  pressure  load  applied  on  the  fluid  domain.  The  time  history  of  the 
charge  is  defined  according  to  Figure  20,  and  it  is  based  on  a  simplified  empirical  model 
[12]. 


Figure  17.  Barge  for  which  vibration  measurements  are  available  [11] 
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Figure  19.  FEA  ABAQUS  model  for  the  barge 


Pressure  vs  Time 


Figure  20.  Definition  of  the  load  from  the  charge  based  on  [12] 

A  velocity  measurement  from  an  UNDEX  test  was  provided  for  a  single  location  placed 
on  the  charge  side  and  close  to  the  centerline  [11].  The  configuration  for  which  the  test 
was  performed  targeted  the  collection  of  data  associated  with  the  performance  of  new 
mounts  used  to  attach  the  two  large  masses  to  the  deck  structure.  The  only  available 
velocity  measurement  corresponded  to  a  single  reference  velocity  measurement  taken 
during  the  test.  The  peak  value  for  the  load  and  the  location  of  the  measurement  point 
comprised  the  random  variables  in  the  probabilistic  model  update  and  the  error  estimation 
analyses.  Based  on  a  suggestion  by  Greg  Harris  (NSWC/IHD),  the  difference  between 
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the  shock  spectra  associated  with  the  measured  and  the  computed  time  histories  of  the 
velocity,  comprised  the  correlation  measurable  in  this  case.  Thus  the  objective  function 
(Equation  (20))  was  modified  to  include  the  difference  in  the  frequency  spectra  rather 
than  the  difference  between  the  time  histories.  The  fast  running  model  represented  by 
equation  (19)  was  still  used  to  compute  the  time  histories  of  the  velocity  during  the  model 
update  and  the  error  estimation  computations,  but  every  time  the  velocity  time  history 
was  processed  in  order  to  compute  the  corresponding  shock  spectrum.  The  velocity  time 
history  prescribed  the  base  excitation  for  a  single  dof  system  with  natural  frequency  “f 
and  the  maximum  value  of  the  velocity  induced  on  the  single  dof  system  was  computed. 
Thus,  the  shock  spectrum  of  maximum  induced  velocity  vs.  the  natural  frequency  of  the 
single  dof  system  is  generated. 

Figure  21  presents  the  area  within  which  the  measurement  location  was  considered  to 
vary.  The  measurement  location  was  considered  as  one  of  the  parameters  with 
uncertainty  because  there  was  no  exact  information  about  the  location  of  the 
measurement  point.  In  Figure  21  the  area  within  which  the  measurement  point  is  allowed 
to  vary  is  highlighted.  The  velocity  values  computed  at  the  nodes  of  the  finite  element 
model  were  interpolated  based  on  the  distance  of  the  measurement  point  under 
consideration  from  the  nodes.  Then  the  corresponding  velocity  time  history  was  used  to 
generate  the  shock  spectrum.  The  transverse  and  the  longitudinal  locations  of  the 
measurement  point  comprised  the  random  variables  expressing  the  uncertainty  in  the 
location  of  the  point.  The  probabilistic  model  update  analysis  was  performed  first.  The 
standard  deviations  for  the  three  random  variables  considered  in  the  analysis  are 
summarized  in  Table  4.  Once  the  model  update  was  completed,  the  error  estimation 
analysis  was  performed.  The  run  time  for  the  model  update  process  was  14  cpu  min  on  a 
Pentium  4  PC,  and  the  run  time  for  the  error  estimation  was  6  cpu  min  on  a  Pentium  4 
PC.  All  the  results  are  summarized  in  Figure  22. 

Uncertainty  in 

Measurement 


Figure  21 .  Uncertainty  in  the  measurement  location 


Random  variables 

Sigma  from  updated  values 

Peak  value  of  load 

1284.0207  psi 

longitudinal  location 

3.0  inch 

transverse  location 

4.5  inch 

Table  4.  Summary  of  the  standard  deviations  for  the  random  variables  considered  in  the 

barge  analysis 
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Shock  spectrum 


The  shock  spectrum  for  the  updated  model  matches  much  better  the  test  data  than  the 
original  numerical  results.  The  boundaries  associated  with  +  or  -  one  standard  deviation 
are  also  presented  in  Figure  22.  It  can  be  observed  that  the  response  spectrum  from  the 
original  model  for  certain  frequencies  is  outside  the  error  bounds  associated  with  one 
standard  deviation  in  the  results.  The  test  data  and  the  numerical  model  fall  within  the  + 
or  -  one  standard  deviation  region  for  all  frequencies. 

This  case  study  demonstrates  that  the  codes  developed  for  probabilistic  model  update 
and  for  error  estimation  can  work  with  real  models  used  in  UNDEX  simulations,  and  that 
they  can  produce  meaningful  results  for  real  applications. 


3.  Potential  developments  during  the  Phase  II  effort 

The  main  objective  for  the  Phase  II  effort  will  be  to  further  develop  the  capabilities  of 
the  probabilistic  model  update  and  the  error  estimation  process  in  order  to  ensure  their 
applicability  in  full  scale  naval  ship  applications.  Further  validation  of  the  new  methods 
should  also  be  pursued.  The  following  tasks  can  be  performed:  (this  list  has  been 
developed  through  interactions  with  NN  and  NSWC/IHD) 

•  Investigate  if  the  Principal  Components  can  be  extracted  from  analyses  other  than 
a  complete  UNDEX  analysis.  This  capability  will  expedite  the  probabilistic 
model  update  process  and  the  error  estimation  when  using  full  scale  ship  models. 
It  will  allow  developing  the  fast  running  models  without  performing  time 
consuming  UNDEX  simulations  but  rather  less  time  consuming  analyses  instead. 

•  Develop  a  variable  screening  process  which  will  permit  to  limit  the  required 
number  of  simulations  at  sample  points  when  developing  the  fast  running  models. 
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•  Establish  a  correlation  metric  between  the  test  data  and  the  numerical  results 
based  on  whether  the  test  data  can  be  recognized  as  a  member  of  the  family  of  the 
numerical  solutions  when  considering  the  uncertainties  in  the  numerical  model. 

•  Apply  the  probabilistic  model  update  and  the  error  estimation  processes  to  an 
automotive  crash  simulation  for  which  both  test  data  and  numerical  models  are 
available. 

•  Work  with  NSWC/IHD  in  order  to  utilize  the  new  processes  developed  by  this 
project  during  the  correlation  phase  between  the  numerical  model  and  the  test  data 
of  their  joint  project  with  the  German  Navy.  A  letter  from  NSWC/IHD  stating 
this  collaboration  with  MES  under  a  Phase  II  contract,  and  endorsing  the  work 
and  the  proposal  from  MES  is  attached. 

•  Work  on  HMMWV  applications  where  the  vehicle  structure  is  subjected  to  the 
blast  load  from  a  mine  or  an  explosive  device. 

•  Develop  an  interface  of  the  probabilistic  model  update  and  error  estimation 
processes  and  the  LMS/Virtual  Lab  program.  This  development  will 
accommodate  commercialization  of  this  technology  beyond  the  Naval 
applications. 
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(Code  334},  Balls ton  Center  Tower  One,  900  Worth  Quincy 
Street,  Arlington,  VA  22217 

Sub j  :  ENDORSEMENT  OF  PHASE  IT  58IR  FOR  PROBABILISTIC  ERROR 

ESTIMATION  IN  MODEL  BASED  PREDICTIONS 


1.  The  Naval  Surface  /Wajpfartt  'Center . 
(N3WC/.IHD)  works  on  the  Navy's  mission-critical  areas  such  as;  . 
energetics  research,-  weapons  product  dbyei  apraerrt  ,>■ . detonation  7 
science r  underwater  warheads,  ordnabce  test ’ and  evaluation,  and 
weapon  simulations.  Specifically  in  the  area  of  accurate 
modelling  and  simulation  of  underwater  explosion  damage  to  naval 
structures,  IHD/HSWC  has  .developed  the  DYS MAS  code.  This 
development  is  based  on  coupling  M  a . hydrocode  with  a 
simulation  model  of  the  target  structure  within  a  parallel 
computing  environment..  In  order  to  validate, and  .demonstrate  th< 
Simulation  capabilities  of  the  OYSMAS  code  IfiD/NSKC  is  currently 
performing  jointly  with  the  German.  Navy  a : full-scale  simulation 
and  test  program  for  a  surface  ship.  :  The . scope  of  this  work  is 
to  model  and  test  eight  different  conflgurafipns  and  compare 
•results  between  simulation  and  test, 

2.  During,  the  past  five  months  Hick  Viahopdulos : of  Michigan 
Engineering  Services,  LLC  .  {MSS L  has-  interacted  with  IHD/NSfcC  anc 
hzs  communicated  results  of  the  work  that  MES  has  performed 
under  the  Phase  I  SBIR  contract . N00014-04-M-0124  titled 
"Probabilistic  error  estimation,  in  model  based  predictions." 

The  technology  that  MES  i.s  developing  is  practical,  technically 
sound,  and  promising  for  combining  results  from  existing  or* 
commercial  non-linear  finite  element  codes  with  information  fret? 
test  in  order  to  update  the  numerical  models  under  uncertainty 
and  provide  an  error  estimation  capability  along  with  the 
numerical  predictions.  The  value  of  their  -  developments  has  beer 
demonstrated  through  case  studies  which  utilize  large: size 
finite  element  models  and  test  data’  from-  Naval  applications. 
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